# -*- coding: utf-8 -*-
from scipy import *
from scipy import optimize

'''
A module to simulate SW particles. The idea is as follow:
First check if the system lies within the Stoner–Wohlfarth astroid
http://en.wikipedia.org/wiki/Stoner%E2%80%93Wohlfarth_astroid
If the system is out of the astroid, there is just one locam minimum. In this case we call the particle saturated. In case of saturation, the particle orientation will be always in the local minimum. If the particle is out of the astroid, we have two local minimum.
The model used in the simulation is as follow:
Totoal energy of the system:
Etot = -A*cos(phi)-cos(theta-phi)^2
Where A=µ0*M*H/K, theta is angle between H and Ku, and phi is between H and M

First the simulation determines the system without thermal effect. 
1. In this case when the system is bistable, the system will be at the state where it follows the slope from the last step.
2. If the system is not bistable, then it just stays at the global minimum.

Then the system will consider the thermal effect. Here the thermal effect f*kT (f=25-30) is compared with the energy barrier:
1. If the system is bistable, the local barrier height Eb has to be considered. if Eb>f*kT, the system simply stays static. Here the effect of thermal fluctuation will in principle change the average value of phi, however, in the limit of Eb>f*kT, only the state that is very close to the local minimum is polulated, and in such case the DOS can be taken as symmetric to a large extend.
2. If the system is bistable and Eb<f*kT, the system will shift to the other local minimum. Then another new barrier Eb2 has to be considered. If Eb2>f*kT, the system stays statically in that state. If Eb2<f*kT, the system is superparamagnetic, and thermal averaging has to be performed.
3. If the system is not bistable, the total barrier height Eb should be compared also with f*kT. and if Eb2>f*kT, thermal averaging has to be done.

The formular for thermal averaging is as follow:

phi_ave = Integrate[Exp[-Etot/kT]*phi*Sin[phi], {phi,0,Pi}]/Integrate[Exp[-Etot/kT]*Sin[phi], {phi,0,Pi}]


To calculate the local minimum, barrier height and the averaging integral, a interplation table is used.
'''

#from func_interp import np1, np2, ne1, ne2, ne3, stable_check

from func_interp2 import np1, np2, ne1, ne2, ne3
from import_math import stable_check


#TESTLIST=[]

#random seed for the monte carlo simulation, so that the result is somehow
#reproducible.
RND_SEED = 1983

# The Boltzmann constant in SI unit
KB = 1.38e-23

#def stable_check(A, theta):
    ## Check for bistable state for SW particle
    ## A hack to revert the theta value into the 0 to pi/2 range.
    #theta = abs(pi*floor(theta/(pi/2))-theta)
    #t = tan(theta)**(1./3.)
    #t2=t*t
    #t4=t2*t2
    #front = sqrt(1+t2+t4)/(1+t2)
    #return A<front

def Etot(A, theta, phi):
    # calculate the Etot
    Etot = -A*cos(phi)-cos(theta-phi)**2
    return Etot
    
def Etotp(A, theta, phi):
    # calculate the Etotp
    Etotp = -sin(2*(theta-phi))+A*sin(phi)
    return Etotp

#int_filename='SW.dat'
## we start to load the data directly
#int_data = loadtxt(int_filename)
#alist = int_data[:,0]
#tlist = int_data[:,1]
#p1list = int_data[:,2]
#p2list = int_data[:,3]
#e1list = int_data[:,4]
#e2list = int_data[:,5]
#e3list = int_data[:,6]

##test_alist = arange(0,0.99,0.04)
##test_tlist = arange(0,pi/2+0.00001,pi/40)
##test_p1list = p1list.copy()
##test_p1list.resize((len(test_alist),len(test_tlist)))
##test_p2list = p2list.copy()
##test_p2list.resize((len(test_alist),len(test_tlist)))

#coord = stable_check(alist,tlist)
#select_alist=alist[coord]
#select_tlist=tlist[coord]
#select_p2list=p2list[coord]
#select_e2list=e2list[coord]
#select_e3list=e3list[coord]


##p1=interpolate.RectBivariateSpline(test_alist,test_tlist,test_p1list)
##p2=interpolate.RectBivariateSpline(test_alist,test_tlist,test_p2list)
#p1=interpolate.interp2d(alist,tlist,p1list,kind='cubic')
#p2=interpolate.interp2d(select_alist,select_tlist,select_p2list,kind='cubic')
#e1=interpolate.interp2d(alist,tlist,e1list,kind='cubic')
#e2=interpolate.interp2d(select_alist,select_tlist,select_e2list,kind='cubic')
#e3=interpolate.interp2d(select_alist,select_tlist,select_e3list,kind='cubic')

## extend the interpolation function using symmetry
#def np1(a, theta):
    #if theta > pi/2:
        #theta = pi - theta
        #out = 2*pi - p1(a, theta)[0]
    #else:
        #out = p1(a, theta)[0]
    #return out

#def np2(a, theta):
    #if theta > pi/2:
        #theta = pi - theta
        #out = 2*pi - p2(a, theta)[0]
    #else:
        #out = p2(a, theta)[0]
    #return out

#def ne1(a, theta):
    #if theta > pi/2:
        #theta = pi - theta
        #out = e1(a, theta)[0]
    #else:
        #out = e1(a, theta)[0]
    #return out

#def ne2(a, theta):
    #if theta > pi/2:
        #theta = pi - theta
        #out = e2(a, theta)[0]
    #else:
        #out = e2(a, theta)[0]
    #return out

#def ne3(a, theta):
    #if theta > pi/2:
        #theta = pi - theta
        #out = e3(a, theta)[0]
    #else:
        #out = e3(a, theta)[0]
    #return out

## sanity check
#x1=0
##x2=0
#for i, a in enumerate(alist):
    #x1 += (p1list[i] - p1(a, tlist[i]))**2
    ##x2 += (p2list[i] - p2(a, tlist[i]))**2
#print x1
##print x2

class SWNPs(object):
    '''A book keepign object to perform simulation based on SW model NPs'''
    def __init__(self, n, theta, phi, volume, Ka, field, magnetisation, temperature, thermal_factor=30):
        # initialize the random number generator
        random.seed(RND_SEED)
        def get(expr, n):
        # to take the possible distribution of theta, phi, volumn and Ka.
            if callable(expr):
                return array(expr(n))
            else:
                return ones(n)*expr
        NP_dict = {}
        NP_dict['theta'] = get(theta, n)
        NP_dict['phi'] = get(phi, n)
        NP_dict['phi2'] = NP_dict['phi'].copy()
        NP_dict['e1'] = zeros(n)
        NP_dict['e2'] = zeros(n)
        NP_dict['e3'] = zeros(n)
        NP_dict['volume'] = get(volume, n)
        NP_dict['Ka'] = get(Ka, n)
        NP_dict['moment_c'] = zeros(n)
        #for i in len(n):
            #NP_list.append(array([get(theta), get(phi), get(volumn), get(Ka)]))
        #NP_list = array(NP_list) #no array at this moment. Not good
        self.NPs = NP_dict # this dict keep track of all the particle information
        self.field = field
        self.magnetisation = magnetisation
        self.temperature = temperature
        self.thermal_factor = thermal_factor
        # calculate the total moment, acts later as the normalization
        totalvolume = NP_dict['volume'].sum()
        totalmoment = totalvolume*magnetisation
        self.totalmoment = totalmoment
        self.relax()
        
    def relaxField(self):
        alist = self.NPs['volume']*self.magnetisation*self.field/(self.NPs['volume']*self.NPs['Ka'])
        bistable = stable_check(alist, self.NPs['theta'])
        self.NPs['bistable']=bistable
        #self.NPs['alist'] = alist
        for i, bis in enumerate(bistable):
            theta = self.NPs['theta'][i]
            if bis:
                pos = self.NPs['phi'][i]
                state1 = (np1(alist[i], theta), ne1(alist[i], theta))
                state2 = (np2(alist[i], theta), ne2(alist[i], theta))
                barrier = ne3(alist[i], theta)
                grad = Etotp(alist[i], theta, pos)
                if state1[0] > state2[0]:
                    highstate = state1
                    lowstate = state2
                else:
                    highstate = state2
                    lowstate = state1
                if (grad<0 and (pos<lowstate[0] or pos>=highstate[0])) or (grad>=0 and pos>=lowstate[0] and pos<highstate[0]):
                    self.NPs['phi'][i] = lowstate[0]
                    self.NPs['e1'][i] = lowstate[1]
                    self.NPs['phi2'][i] = highstate[0]
                    self.NPs['e2'][i] = highstate[1]
                    self.NPs['e3'][i] = barrier
                else:
                    self.NPs['phi'][i] = highstate[0]
                    self.NPs['e1'][i] = highstate[1]
                    self.NPs['phi2'][i] = lowstate[0]
                    self.NPs['e2'][i] = lowstate[1]
                    self.NPs['e3'][i] = barrier
                #TESTLIST.append(theta)
                #pos1 = np1(alist[i], theta)
                #pos2 = np2(alist[i], theta)
                ##pos1_test = array([pos1-2*pi, pos1, pos1+2*pi])
                ##pos2_test = array([pos2-2*pi, pos2, pos2+2*pi])
                ##if theta > pi/2: 
                ##    print 'theta = %f' %theta
                ##    print 'pos2 = %f' %pos2
                #distance = pos1 - self.NPs['phi'][i]
                #if abs(distance) > pi:
                    #distance = -sign(distance)*(2*pi-abs(distance))
                #test=distance*Etotp(alist[i], self.NPs['theta'][i],self.NPs['phi'][i])
                ## The test is as follow. The phi should not go pass any energy barrier, 
                ## thus the direction of movement always decrease the energy. This means:
                ## If phi increase, energy decrease, or Etotp negative
                ## If phi decrease, energy decrease, Etotp positive.
                ## Thus test should always be negative
                #if test<=0:
                    #self.NPs['phi'][i] = pos1
                    #self.NPs['e1'][i] = ne1(alist[i], theta)
                    #self.NPs['phi2'][i] = pos2
                    #self.NPs['e2'][i] = ne2(alist[i], theta)
                    #self.NPs['e3'][i] = ne3(alist[i], theta)
                #else:
                    #self.NPs['phi'][i] = pos2
                    #self.NPs['e1'][i] = ne2(alist[i], theta)
                    #self.NPs['phi2'][i] = pos1
                    #self.NPs['e2'][i] = ne1(alist[i], theta)
                    #self.NPs['e3'][i] = ne3(alist[i], theta)
            else:
                pos = np1(alist[i], theta)
                self.NPs['phi'][i] = pos
                self.NPs['phi2'][i] = pos

    def relaxTemp(self):
        thermal_energy = self.temperature*KB
        for i,bis in enumerate(self.NPs['bistable']):
            if bis:
                barrier_energy = self.NPs['Ka'][i]*self.NPs['volume'][i]*(self.NPs['e3'][i] - self.NPs['e1'][i])
                if self.thermal_factor*thermal_energy > barrier_energy:
                    # if the particle is bistable, but with enough thermal energy, it will reach
                    # thermal equilibrium between two states.
                    teste1 = self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e1'][i]
                    teste2 = self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e2'][i]
                    if teste1 <= teste2:
                        factor = exp(-(teste2-teste1)/thermal_energy)
                        moment_contribution = self.magnetisation*self.NPs['volume'][i]*(cos(self.NPs['phi'][i])+factor*cos(self.NPs['phi2'][i]))/(1+factor)
                    else:
                        # temperature induced flip
                        t = self.NPs['phi'][i]
                        self.NPs['phi'][i] = self.NPs['phi2'][i]
                        self.NPs['phi2'][i] = t
                        t = self.NPs['e1'][i] 
                        self.NPs['e1'][i] = self.NPs['e2'][i]
                        self.NPs['e2'][i] = t
                        barrier_energy = self.NPs['Ka'][i]*self.NPs['volume'][i]*(self.NPs['e3'][i] - self.NPs['e1'][i])
                        if self.thermal_factor*thermal_energy > barrier_energy:
                            teste1 = self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e1'][i]
                            teste2 = self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e2'][i]
                            factor = exp(-(teste2-teste1)/thermal_energy)
                            moment_contribution = self.magnetisation*self.NPs['volume'][i]*(cos(self.NPs['phi'][i])+factor*cos(self.NPs['phi2'][i]))/(1+factor)
                        else:
                            moment_contribution = self.magnetisation*self.NPs['volume'][i]*cos(self.NPs['phi'][i])
                    #factor1=exp(-self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e1'][i]/thermal_energy)
                    #factor2=exp(-self.NPs['Ka'][i]*self.NPs['volume'][i]*self.NPs['e2'][i]/thermal_energy)
                    #moment_contribution = self.magnetisation*self.NPs['volume'][i]*(factor1*cos(self.NPs['phi'][i])+factor2*cos(self.NPs['phi2'][i]))/(factor1+factor2)
                    ## temperature induced flip
                    #if factor2 > factor1:
                        ## exchange the phi and e
                        #t = self.NPs['phi'][i]
                        #self.NPs['phi'][i] = self.NPs['phi2'][i]
                        #self.NPs['phi2'][i] = t
                        #t = self.NPs['e1'][i] 
                        #self.NPs['e1'][i] = self.NPs['e2'][i]
                        #self.NPs['e2'][i] = t
                else:
                    moment_contribution = self.magnetisation*self.NPs['volume'][i]*cos(self.NPs['phi'][i])
                self.NPs['moment_c'][i] = moment_contribution
            else:
                moment_contribution = self.magnetisation*self.NPs['volume'][i]*cos(self.NPs['phi'][i])
                self.NPs['moment_c'][i] = moment_contribution
                
    def relax(self):
        self.relaxField()
        self.relaxTemp()
        
    def getMoment(self):
        m_tot = sum(self.NPs['moment_c'])
        m_tot = m_tot / self.totalmoment # normalize to saturation.
        return m_tot
        
    def updateField(self, field):
        self.field = field
        self.relax()
    
    def updateTemperature(self, temperature):
        self.temperature = temperature
        self.relaxTemp()
        

# Test the object.
# All values in SI unit.
#testvolume = (9e-9)**3
#testKa = 1e4
#testphi = array([0, pi/2, pi, 3*pi/2])
#testtheta = array([0, pi/6, pi/4, pi/3])
#testn = 4
#testmag = 1700e3
#testfield = 0.01

#testSW = SWNPs(testn, testtheta, testphi, testvolume, testKa, field=testfield, magnetisation=testmag, temperature = 100)

#print testSW.getMoment()


###
# now define functions to generate random initial conditions.
def generate_theta(n=None):
    '''generate random number in 0 to pi range, with sin(theta) weighting.'''
    if n:
        pool = random.random(n)
    else:
        pool = random.random()
    out = arccos(1-2*pool)
    return out

def generate_phi_old(n=None):
    '''generate phi in 0 to pi, with sin(phi) weighting.
    In order to keep the numerical error in track, we make a hack, that the phi has equal weigh
    '''
    if n:
        pool = random.random(n)
        out = arccos(1-2*pool)
    else:
        pool = random.random()
        out = arccos(1-2*pool)
    return out

def generate_phi(n=None):
    '''generate phi in 0 to 2pi, with sin(phi) weighting.
    In order to keep the numerical error in track, we make a hack, that the phi has equal weigh
    
    Note that the n has to be a odd number, in order this works.
    '''
    if n:
        m = int(n/2)
        pool = random.random(m)
        down = arccos(1-2*pool)
        up = down + pi
        out = concatenate((up, down))
    else:
        pool = random.random()
        if pool >=0.5:
            out = 2*pi-arccos(4*pool-3)
        else:
            out = arccos(1-4*pool)
    return out

def generate_volume(n=None):
    '''generate the volume distribution for a spherical particle with log-normal dist.
    Not the average and variance of the log-normal distribution is NOT the mu and sigma factor here.
    One should calculate the value again.
    
    Use the same hack as the above function, in order to reduce numeric error
    '''
    mu = 2.58041
    sigma = 0.121255
    m = n/2
    dia = random.lognormal(mu,sigma,m)*1e-9
    volume = pi*dia**3/6
    volume = concatenate((volume, volume))
    return volume

def simulate_zfc_fc(NP, field, temp_list):
    '''The function to perform a ZFC-FC simulation. The NP should first be initialized.
    The function performs the operation on the NP object.
    The result of the simulation will be generated.
    '''
    temp_list.sort()
    NP.updateField(field)
    ZFC_moment = []
    for temp in temp_list:
        NP.updateTemperature(temp)
        ZFC_moment.append(NP.getMoment())
    ZFC_moment = array(ZFC_moment)
    FC_moment = []
    for temp in temp_list:
        NP.updateTemperature(temp)
        FC_moment.append(NP.getMoment())
    FC_moment = array(FC_moment)
    return ZFC_moment, FC_moment

def fit_zfc_fc_ka(ka, templist, ZFCdata, FCdata, N, theta, phi, volume, magnetisation, field):
    # The ZFC and FC data has to be normalized. Else it will NEVER work!!!!
    # The data is normalized so that the max of all data is 1, and min of all is 0
    FitNP = SWNPs(N, theta, phi, volume, ka, 0, magnetisation, 0.00001)
    templist.sort()
    ZFCfit, FCfit = simulate_zfc_fc(FitNP, field, templist)
    norm_up = max(max(FCfit), max(ZFCfit))
    norm_down = min(min(FCfit), min(ZFCfit))
    #ZFCfit = (ZFCfit-norm_down)/(norm_up-norm_down)
    #FCfit = (FCfit-norm_down)/(norm_up-norm_down)
    ZFCfit = ZFCfit/norm_up
    FCfit = FCfit/norm_up
    err = concatenate((ZFCdata-ZFCfit, FCdata-FCfit))
    print 'estimate error is: %f' %(err**2).sum()
    return err

def fit_zfc_fc_ka_mag(x0 , templist, ZFCdata, FCdata, N, theta, phi, volume, field):
    # The ZFC and FC data has to be normalized. Else it will NEVER work!!!!
    # The data is normalized so that the max of all data is 1, and min of all is 0
    ka, magnetisation = x0
    FitNP = SWNPs(N, theta, phi, volume, ka, 0, magnetisation, 0.00001)
    templist.sort()
    ZFCfit, FCfit = simulate_zfc_fc(FitNP, field, templist)
    norm_up = max(max(FCfit), max(ZFCfit))
    norm_down = min(min(FCfit), min(ZFCfit))
    ZFCfit = (ZFCfit-norm_down)/(norm_up-norm_down)
    FCfit = (FCfit-norm_down)/(norm_up-norm_down)
    #ZFCfit = ZFCfit/norm_up
    #FCfit = FCfit/norm_up
    err = concatenate((ZFCdata-ZFCfit, FCdata-FCfit))
    print 'use ka=%f ms=%f estimate error is: %f' %(ka,magnetisation,(err**2).sum())
    return err

def fit_zfc_fc_ka_size(x0 , templist, ZFCdata, FCdata, N, theta, phi, magnetisation, field):
    # fit the ZFCFC with Ka and size distribution
    ka, mu, sigma = x0
    def generate_lognorm(n=None):
        '''generate the volume distribution for a spherical particle with log-normal dist.
        Not the average and variance of the log-normal distribution is NOT the mu and sigma factor here.
        One should calculate the value again.
        
        Use the same hack as the above function, in order to reduce numeric error
        '''
        #mu = 2.58041
        #sigma = 0.121255
        #global mu
        #global sigma
        m = n/2
        dia = random.lognormal(mu,sigma,m)*1e-9
        volume = pi*dia**3/6
        volume = concatenate((volume, volume))
        return volume
    
    FitNP = SWNPs(N, theta, phi, generate_lognorm, ka, 0, magnetisation, 0.00001)
    templist.sort()
    ZFCfit, FCfit = simulate_zfc_fc(FitNP, field, templist)
    norm_up = max(max(FCfit), max(ZFCfit))
    norm_down = min(min(FCfit), min(ZFCfit))
    ZFCfit = (ZFCfit-norm_down)/(norm_up-norm_down)
    FCfit = (FCfit-norm_down)/(norm_up-norm_down)
    #ZFCfit = ZFCfit/norm_up
    #FCfit = FCfit/norm_up
    err = concatenate((ZFCdata-ZFCfit, FCdata-FCfit))
    print 'use ka=%f mu=%f estimate error is: %f' %(ka,mu,(err**2).sum())
    return err


def fit_zfc_fc_ka_size_para(x0 , templist, ZFCdata, FCdata, N, theta, phi, magnetisation, field):
    # fit the ZFCFC with Ka and size distribution
    ka, mu, sigma, para = x0
    if sigma < 0: sigma = 1e-5
    if mu < 0: mu = 1e-5
    def generate_lognorm(n=None):
        '''generate the volume distribution for a spherical particle with log-normal dist.
        Not the average and variance of the log-normal distribution is NOT the mu and sigma factor here.
        One should calculate the value again.
        
        Use the same hack as the above function, in order to reduce numeric error
        '''
        #mu = 2.58041
        #sigma = 0.121255
        #global mu
        #global sigma
        m = n/2
        dia = random.lognormal(mu,sigma,m)*1e-9
        volume = pi*dia**3/6
        volume = concatenate((volume, volume))
        return volume
    
    FitNP = SWNPs(N, theta, phi, generate_lognorm, ka, 0, magnetisation, 0.00001)
    templist.sort()
    ZFCfit, FCfit = simulate_zfc_fc(FitNP, field, templist)
    ZFCfit = ZFCfit + para/templist
    FCfit = FCfit + para/templist
    norm_up = max(max(FCfit), max(ZFCfit))
    norm_down = min(min(FCfit), min(ZFCfit))
    ZFCfit = (ZFCfit-norm_down)/(norm_up-norm_down)
    FCfit = (FCfit-norm_down)/(norm_up-norm_down)
    #ZFCfit = ZFCfit/norm_up
    #FCfit = FCfit/norm_up
    err = concatenate((ZFCdata-ZFCfit, FCdata-FCfit))
    print 'use ka=%f mu=%f simga=%f para=%f estimate error is: %f' %(ka,mu,sigma,para,(err**2).sum())
    return err



if __name__ == '__main__':
    # For test only
    #TestNP = SWNPs(5000, generate_theta, generate_phi, generate_volume, 1e4, 0, 1700e3, 0.00001)
    #temp_list = arange(5,300,5)
    #zfc, fc = simulate_zfc_fc(TestNP, 1e-3, temp_list)
    #out = array([temp_list, zfc, fc]).transpose()
    #savetxt('testzfcfc', out)
    
    zfcfile = 'ZFC-set1.dat'
    fcfile = 'FC-set1.dat'
    #Ka_start = 3e4
    #Ka_start = 5.4e4
    Ka_start = 40000
    #mu_start = 2.58041
    #sigma_start = 0.121255
    N = 10000
    #magnetisation = 1700e3
    #magnetisation = 169815
    magnetisation = 170000
    #para_start = 5e-16
    field = 0.01
    ZFCdatas = loadtxt(zfcfile)
    FCdatas = loadtxt(fcfile)
    templist = ZFCdatas[:,0]
    ZFCdata = ZFCdatas[:,1]
    FCdata = FCdatas[:,1]
    # normalize the data
    norm_up = max(max(FCdata), max(ZFCdata))
    norm_down = min(min(FCdata), min(ZFCdata))
    ZFCdata = (ZFCdata-norm_down)/(norm_up-norm_down)
    FCdata = (FCdata-norm_down)/(norm_up-norm_down)
    #ZFCdata = ZFCdata / norm_up
    #FCdata = FCdata / norm_up
    #start the simualtion
    #res = optimize.leastsq(fit_zfc_fc_ka, Ka_start,args=(templist, ZFCdata, FCdata, N, generate_theta, generate_phi, generate_volume, magnetisation, field))
    res = optimize.leastsq(fit_zfc_fc_ka_mag,array([Ka_start, magnetisation]),args=(templist, ZFCdata, FCdata, N, generate_theta, generate_phi, generate_volume, field))
    #res = optimize.leastsq(fit_zfc_fc_ka_size,array([Ka_start, mu_start, sigma_start]),args=(templist, ZFCdata, FCdata, N, generate_theta, generate_phi, magnetisation, field))
    #res = optimize.leastsq(fit_zfc_fc_ka_size_para,array([Ka_start, mu_start, sigma_start, para_start]),args=(templist, ZFCdata, FCdata, N, generate_theta, generate_phi, magnetisation, field))
    
    print res
    Ka = res[0][0]
    magnetisation = res[0][1]
    #mu = res[0][1]
    #sigma = res[0][2]
    #para = res[0][3]
    
    def generate_lognorm(n=None):
        '''generate the volume distribution for a spherical particle with log-normal dist.
        Not the average and variance of the log-normal distribution is NOT the mu and sigma factor here.
        One should calculate the value again.
        
        Use the same hack as the above function, in order to reduce numeric error
        '''
        #mu = 2.58041
        #sigma = 0.121255
        #global mu
        #global sigma
        m = n/2
        #if sigma < 0: sigma = 1e-5
        #if mu < 0: mu = 1e-5
        dia = random.lognormal(mu,sigma,m)*1e-9
        volume = pi*dia**3/6
        volume = concatenate((volume, volume))
        return volume
    
    FitNP = SWNPs(N, generate_theta, generate_phi, generate_volume, Ka, 0, magnetisation, 0.00001)
    ZFCfit, FCfit = simulate_zfc_fc(FitNP, field, templist)
    #ZFCfit = ZFCfit + para/templist
    #FCfit = FCfit + para/templist
    ZFCfit = array([templist,ZFCfit]).transpose()
    FCfit = array([templist,FCfit]).transpose()
    savetxt('ZFCfit',ZFCfit)
    savetxt('FCfit',FCfit)

    
